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Abstract 

In this paper, we study the mechanics of statistically non-uniform two-phase elastic 
discrete structures. In particular, following the methodology proposed in (Luciano 
and Willis, Journal of the Mechanics and Physics of Solids 53, 1505-1522, 2005), en- 
ergetic bounds and estimates of the Hashin-Shtrikman- Willis type are developed for 
discrete systems with a heterogeneity distribution quantified by second-order spatial 
statistics. As illustrated by three numerical case studies, the resulting expressions 
for the ensemble average of the potential energy are fully explicit, computationally 
feasible and free of adjustable parameters. Moreover, the comparison with reference 
Monte-Carlo simulations confirms a notable improvement in accuracy with respect 
to approaches based solely on the first-order statistics. 

Key words: inhomogeneous material (B); structures (B); energy methods (C); 
probability and statistics (C) 



1. Introduction 

Discrete material models, which represent a material as a network of particles 
interacting via inter-particle potentials, have received a steadily increasing atten- 
tion in the fields of theoretical, co mputational and applied materials sci e nce i n 



the last dec a de, see , e.g., reviews by lOstoja-Starzewskil (120021 ): lAlava et al.l (120061 ) 



Blanc et al.l (j2007aj) and references therein. From the engineering point of view. 
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the interest has been nourished by the possibihty to address, in a conceptually 
simple framework, the interplay among the intrinsic material heterogeneities, dis- 
creteness and randomness on different levels of resolution. Examples of the suc- 
cessful application of discrete mode ls include the simulation of materials with com- 



plex microstructures such as paper (lOstoia-Starzewski and Stahll . 1200 It iBronkhorst 



2003 ). biological materials ( Arnoux et al.l. 120021) , low-d ensity materials (jChristensen 



2000) and granular media (IMiehe and Dettmarl . l2004l ) . Another field of application 



where the discrete modeling concept plays an irreplaceable role is the analysis of 
localized phenomena in het erogeneous media, such as local buckling and delam- 



ination in thin films, e.g. ( jjaglal . 120071 : IVellinga et al.l . 120081 ) . or, most typically, 



the simulation of damage and fracture processes in cohesive-fricti o nal materials 



see ( 


Lilliu and van Mier 




2003: 


Ibrahimbeeovic and Delaulace 


2006 


: Chudoba et al.. 12006: 


Vofechovskv and Chudoba. 


2006h 



The closely related problem of establishing a rigorous link between a discrete rep- 
resentation and its equivalent continuum response has been the focus of numerous 
works. The goal of such studies is to identify an appropriate continuum represen- 
tation with materials constants directly related to the underlying discrete system. 
Within the computational approaches, perhaps the most prorn i nent example is the 
'local' Quasi- Continuum method, introduced by iTadmor et al.l (119961 ). in which the 
effective behavior of a material point is adaptively deduced fro m the response of its 



Miller and Tadmor 



finite neighborhood, constrained by the Cauchy-Born rule; see 
(I2OO2I ) for further details and discussion on related concepts. 

Complementary to the numerical treatment, a notable research effort has been put 
into a rigorous interpretation of infinite-size limits of discrete models from the point 
of view of standard and generalized continua. For systems interacting via potentials 
satisfying suitable growth conditions, a general local cont inuum representation is 
currently available, proven using the tools of F-converge nce (lAlicandro and Cicalesd . 



20051 ) or the thermodynamical limit procedure due to iBlanc et al.l (|2002[ ). These 



results were further utilized to pr ovide rigorous bounds 1 1 the effective macroscopic 
conductivity of discrete lattices ( iBraides and Francfortl . l2004j ) or as a theoretical 
support fo r the quasi- con t inuum approximation to ab initio calculations of material 



constants ( iGavinia et al.l . 120071 ). Re cently, both frameworks wer e suc cessfully ex- 



tended to the stochastic setting, see lAlicandro et al.l (120071 l2010l ) and iBlanc et al. 
(l2007bl H). In addition, the y a lidity of t he Cauchy-Born rule wer e rigor ously ex- 
amined in iFriesecke and Theill ( 120021 ) and iBerezhnyy and Berlyandl (120061 ) for both 
regular and irregular networks, thus explicitly demonstrating potential limitations of 
the Cauchy-Born type continuum approximation when applied to discrete localized 
phenomena. 
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In such cases, generalized continuum theories provide a well-estabhshed way to 
introduce an internal lengthscale to the problem, thereby preventi ng the pathological 
local i zation or singularities of rnechanical fields, see e.g. reviews (IGanghoffer et al.l . 
19991 : iBazant and Jirasekl . 120021 : lEringenl . |2002| ). Particular examples of gradient- 



based theories include the one- and two-dim ensi onal large strain elasticity stud- 
i es by iTriantafyllidis and BardenhagenI ( 1l993l ) and iBardenhagen and TriantafvUidis 
199J), the micropolar continuum description of lPradel and Sa bl fll998h and' 
fl2007l) or the arbitrary-order convex expansion scheme due to lArndt and Griebel 



Martinsson and Babuska 



fl2005h . Finally, an exhaustive analysis of one-dimensional systems with generic 
nearest-neighbor interactions rigorously demonstrated that the lim it behavior may 
exhib it both diffuse as well as localized crack ing, in the deterministic ( iBraides and Gelli 
2002) and stochastic (jlosifescu et al.l . 120011) setting, includi ng the numerical a nalysis 



of discrete-to-continuum coupling in the deterministic case flBlaiic 



recently, these results were extended by lBraides and Piatnitskil ( l2008l ). who studied 



eta 



20051). More 



continuum limits of lattices with randomly distributed defects and showed that the 
effective behavior is governed by percolation phenomena. The treatment of finite- 
size discrete systems is, to our best knowledge, much less developed and is typically 
limited to f itting of phenornenological constitutive relations to numerical simulation 
results, cf. (IRinaldi and Lail . 120071 : iGrassl and Jirasekl . l2010l . and references therein). 

In the present paper, we address in detail a specific problem of the mechanics of 
random discrete media, namely the formulation of total potential energy estimates for 
finite binary lattices with a fixed geometry and a heterogeneity distribution described 
by second-order spatial statistics. Variational bounds and estimates are established 



by extending the recent works of Luciano and Willis ( 2005 . 2006! ) related to the 



Galerkin d iscretization of the stochastic Hashin-Shtrik man- Willis (HSW) variational 
principles ( iHashin and Shtrikmanl . Il962l : IWillij . Il977l ). Our motivation for focusing 
on finite-sized systems and the potential energy instead of the more common con- 
tinuum setting and local stress- or strain-related quantities arises from the following 
considerations: 



The separation-of-scales assumption is inherently not applicable when dealing 
with finite discrete structures. This renders the resulting theory well-suited to 
predict the statistics of localized responses. 



In view of recent advance s in variationa l models of complete damage (iBouchitte et al. 



20091 : iMielke et al.l . l2010l : iMielkd . 120091 ). the global energetic bounds/estimates 
provide an essential ingredient for the development of 'rational' damage me- 
chanics of discrete networks. 
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Due to the simple structure of the underlying theory, the relevant statistics can 
be characterized with a generality which is currently not available for continu- 
ous systems, cf. Section 12.21 

Highly accurate estimates of the quantities of interest for general structures and 
loading regim es can be determined on the basis of simple Monte-Carlo simula- 



tions, see also ISharif-Khodaei and ZemanI (j2008[ ) for a related one-dimensional 



study in the continuous setting. 

The remainder of the paper is organized as follows. In Section [21 the relevant 
steps of the problem definition are specified for both deterministic and randomized 
systems. Energetic bounds and estimates are derived in Section [3l In Section HI 
results of numerical studies are presented to assess their accuracy and limitations. 
Finally, Section [5] collects concluding remarks and comments on future extensions of 
the method. 



2. Problem description 



This Section is devoted to the problem statement, starting with a brief summary 
of structural statics, in order to introduce our notation, followed by the specification 
of the stochastic framework an d quantities of interest. St andard notations and results 
of linear algebra are employed (IHorn and Johnsonl . ll990l ). with a, a and A denoting a 
scalar quantity, a vector (column matrix) and a generic matrix, respectively. Matrix 
indexing is used when appropriate, i.e. given two index sets i and j with cardinalities 
\i\ and '^^A E RI^I^I^'I denotes the appropriate sub-matrix of A, while gives 
the corresponding matr ix rows and ^'A"*" abbrevia tes Moreover, the matrix 

formalism developed in iJirasek and BazantI (120011 ) for general discrete structures is 
systematically adopted. 



2.1. Summary of discrete media mechanics 

Consider a discrete structure consisting of Nn nodes with coordinates Xi G M'^, 
i = {1, 2, ... , A^„} and d G {2, 3}, connected by A^e discrete elements. On the level 
of a single element e G {1,2,..., N^}, the generalized kinematic equations take the 
form 



Bpdp 



where Be G M^'' is the vector of generalized strains, the vector dg G R^^'' contains the 
Nd generalized displacements at both element nodes and Be G M^^^^^^ denotes the 
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element kinematic matrix. The corresponding generalized element stresses Se G 
then follow from 

Se = DeBe, (2) 

with Z?e G R-^^^^" denoting a matrix of generalized material stiffness. On the 
structural level, the relations ([1]) and (|2]) can be assembled into the form 

e = Bd, s = De, 

where, e.g., B E R^^^-x^^^d^ D e R^-^sxn,Ns ^^d d E R^"^'^ stand for the global 
kinematic matrix, block-diagonal generalized stiffness matrix and displacement vec- 
tor defined as 



B = A Be, D = A De and d = A de, (3) 

e=l e=l e=l 

with the symbol A representing the assembly operation, cf. f Jirasek and Bazantl . 
200 ll . Appendix A). The remaining matrices and vectors are defined analogously. 

In order to specify kinematic constraints on the structure, we partition the prob- 
lem degrees of freedom (DOFs) into two sets 

cUu = {1,2,...,A^„A^,}, cnu = (D, ker(^"B) = {0}, (4) 

where c and u collect the known (constrained) and unknown DOFs and ker(A) 
denotes the kernel of a matrix A. Note that the last condition in Eq. enforces 
the elimination of rigid-body modes. When subjecting the structure to an additional 
nodal load G M'"' acting on free DOFs, the unknown displacements ^d can be 
found by solving the unconstrained quadratic optimization problem 

^d = arg min E{d), 



Ne 

A De 

e=l 



Ne 

A de 

e=l 



where the total potential energy function E : 



u 



Bid) 



d^ ^d"^ 



UCJ^ 



— )■ M is provided by 





' d ' 







(5) 



(6) 



ucj^ = ^ B^D ''B G RI""!^!"! being a sub-matrix of the global 
= B^DB. The symbol ' 'arg min" appearing in Eq. (E]) denotes 



with, for example, 
stiffness matrix K 
the minimizer of the objective function verifying 

E("d) < E{d) Vd G 



(7) 



where the equality is attained only for the test displacement d coinciding with the 
true solution due to positive definiteness of '^'^K. The optimality conditions for '"d 
then yield the global equilibrium equations in the form 

'*"li:"d = "/-^^K'^d. (8) 
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2.2. Stochastic setting 

We now proceed with the introduction of a suitable framework for binary random 
discrete media, i.e. structures in which every element can be found in one of two 
distinct states r G {1,2}. Due to the discrete nature of the problem at hand, the 
ensemble space § collecting all structural configurations is finite-dimensional and as 
such can be enumerated using an index a, 

aG§={l,2,...,2^^}. (9) 

The complete statistical characterization of the discrete stochastic system is then 
simply provided by assigning probabilities to individual configurations a stored 
in the probability distribution vector 

r isi ] 

HEA=<fle Ml^l,/I(a) > Va e = 1 > . (10) 

The ensemble average of a configuration- dependent quantity /(«) for a given prob- 
ability distribution ^ G A is defined as 

|S| 

(/);. = E/(«)M«)- (11) 

a=l 

Of particular importance is the state characteristic vector x^^K^t) defined via 

(r). X _ J 1 if element e is in state r for configuration a, , . 

Xe = I otherwise, ^ ' 

quantifying the spatial distribution of individual states in a given configuration a. 
For the current, binary, case, X^^K^i) can be explicitly expressed in the form 

X^'\a) = {a- 1)b, x'^'\a) + x^'\a) = 1, (13) 

where Um provides the value of a natural number n in the binary notation using 
A'^e digits arranged in a column matrix, see Figure [H for an illustration. 

Following the analogy with the quantifi cation of the spatial statistics of random 
heterogeneous media, e.g. (ITorquatol . |2002| ) . we introduce an A^'e x two-unit prob- 
ability matrix related to a given probability distribution fi in the form 

pirs) ^ (;^M^WT>^ = 5^x(^n«)X^^n«)V(«), (14) 
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where an individual entry PjfJ^ represents the probabihty of states r and s being 
assigned to elements e and e' (note that the explicit dependence on /x is dropped for 
the sake of notational brevity). The two-unit matrices related to distinct states are 
not independent of each other: using f lT^ o yields 

p(12) ^ p(l)lT_p(ll), (15) 
p(21) ^ lp(l)T_p(ll), (16) 
P(22) = llT_p(l)lT_lp(l)T^p{ll)^ (17) 

with p^^^ = diag(P'-''^''). It is therefore sufficient to concentrate on the statistics P^^^^ 
below. 

It directly follows from the definition f[T^ that any two-element probability matrix 
has to be located in a convex hull of rank-one products of the characteristic vectors 
of individual configurations: 

p(ii) eMQ = conv {x^'\l)x^'\l)\ x^'\2)x^'\2)\ . . . , x^^HlSD^^^HlSI)"'} , (18) 

coinci ding with the Boolean quadratic polytope completely characterized by lPadberg 



(119891 ). Conversely, two-element probability matrices can be used as a convenient 
re-parameterization of A. To that end, we introduce a set storing all probability 
distributions compatible with matrix P^^^^ as 

M(p(^i)) = {neA, p(i^) = ix^'V^'')^} , (19) 

thus establishing a partial statistical characterization when full information is not 
availablell] 



^To the best of our knowledge, no such results are available for general multi-unit probability 
functions. In addition, the treatment of higher-order statistics leads to a substantial increase 
of storage and computing requirements. Therefore, we limit our attention to the second-order 
framework and leave its extension to future work. 

a = 1 a = 2 

1 ^kk.O O ^2 O O ^ 1 » O C^^H 2 

X^'Hl) = [0,0]T,;^(2)(i) = [i,i]T ;^(i)(2) = [0, 1]T,;^(2)(2) = [1,0]^ 

a = 3 a = 4 

Z(')(3) = [l,0]"r,;fP)(3) ^ [o^ijT ^(i)(4) ^ [i,i]T^;^(2)(4) ^ [o^ojT 

Figure 1: Ensemble space and state characteristic vectors for a two-element structure; state r = 1 
corresponds to a black element, r = 2 to a gray element. 
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At this point, we can introduce the terminology used hereafter. Given a two- unit 
probabihty matrix P^^^\ the associated discrete system is called statistically uniform 
if 



(1) 



P, 



(11) 



P(x, 



X' 



for e, e = 1, 2, 



(20) 



i.e. the one-unit probabilities of individual elements are independent of e, while the 
two-unit probabilities depend only on the difference of coordinates of element centers 
Xe and ajg. In addition, a function /(P*-^^-*) is statistically local (or first-order) if it 
depends on the first-order statistics only: 



(1)^ 



(21) 



In the opposite case, we use the adjectives statistically non-local (or second-order). 

It is perhaps instructive to briefly comment on similarities and differences with the 
analogous development in the continuous case. In particular, the two-element prob- 
ability matrix is a direct analogue of the two-point probabil ity functions for generic 



statis tically non-uniform and anisotropic microstructures (IWillisl . Il98ll : iTorquatd . 
2OO2I ). The set Bq then coincides with the set of all admissible two-point prob- 
ability functions, the analysis of which has recently received substantial attention 
in the field of random heterogeneous materials, eventually resulting in a corn plete 
char acterization for st atistically uniform media as shown by iQuintanillal (120081 ): see 
also IJiao et al.l ( 120071 ) for further discussion. In this connection, the advantage of 
consideri ng finite d i screte systems becomes immediately visible, since the character- 
ization of iPadbergI (119891 ) is free of assumption of statistical uniformity. 



2.3. Statics of two-phase random lattices 

When assigning specific material properties to each state r, the previously intro- 
duced framework can be readily adopted to two-phase elastic heterogeneous lattices. 
In particular, the configuration-dependent material stiffness matrix on the element 
level, De{a), can be written as 



(r) 
e 5 



(22) 



r=l 



where D^j'' are generalized material stiffness matrices of individual states. On the 
structural level, the stiffness distribution is described by 



^ r=l ^ r=l 



a 



(23) 
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where a • A denotes a block Hadamard-like product implementing the assembly 
operation. 

Consider now the response of a discrete structure with a stochastic configuration- 
dependent generalized material stiffness matrix D{a) subject to deterministic load- 
ing conditions specified in terms of prescribed displacements and generalized nodal 
forces For each configuration a G S, the energy minimizer is defined as 



''d{a) = argjnin E{d;a) 



(24) 



with the stored energy function introduced analogously to ([6]): 



E(d:a) 



d^ ^d"^ 



d 



d^ 



(25) 



The ensemble average of the optimal energy for a given probability distribution fi is 
then simply a weighted sum 



|§| 



{Erd))^ = J2ECdia))f,ia). 



(26) 



a=l 



A full specification of the probability distribution is, however, rarely available 
and even with complete information at hand, evaluating ( 126|) requires the solution 
of 2^" problems, which becomes rapidly unfeasible even for moderate-size problems. 
Therefore, we rely on a partial statistical characterization in terms of two-element 
probabilities and attempt to establish energetic bounds/estimates in the form 

n-{p^^^^) < H^(p(")) ^ {E^d))^ < n+{p^^^^) v/i e M(p(^i)), (27) 

refiecting the limited probabilistic characterization. In addition to the energetics, we 
also provide elementary statistics of the nodal displacements related to the bounds 
and estimates of the energy introduced in Eq. ([2] 



3. Hashin-Shtrikman- Willis-type estimates 

In this Section, we provide explicit energetic bounds for random finite-size net- 
works by reconsidering the Hashin-Shtrikman variational principles for random het- 
erogeneous media in the discrete setting. To make the exposition more readable, the 
derivations are structured in six consecutive steps. 
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3.1. Reference structure and general ized 'polarization stresse s 



Following the conceptual lead of iHashin and ShtrikmanI (jl962[ ). we introduce a 
reference deterministic structure characterized by a positive-definite generalized ma- 
terial stiffness matrix ZJ'-^^ and consider a realization-dependent quadratic form 



I 



T 

e 



(28) 



where e G M^'=^^ is a test generalized strain vector and the auxiliary variable r G 
]^Ne Ns commented o n later. By virtue of the Schur complement lemma, 

cf. ( Horn and Johnsonl . 1990 . Section 7.7.6), the form f l28|) is positive-semidefinite as 
long as [D^^^ — D{a)) is positive definite, leading to a bound 



with the equality reserved for r equal to 

'r>(°) - Dia 



T a 



1^T(^(0) 



s[a] 



Dia)) 



-1; 



(29) 



(30) 



which minimizes the right hand side of f l29|) for a given realization a and test gen- 
eralized strain field e. The variable r therefore corresponds to a generalized po- 
larization stress_assod^ied__w^^ the re ference stiffness matrix Z)*-"-* and generalized 
strain e (IHashin and ShtrikmanI . Il962l ) . 

Inequality (129!) leads, for an arbitrary positive-definite [D^^^ —D{a)), to an upper 
bound on the configuration-dependent internal energy of the structure expressed in 
terms of an auxiliary structure with stiffness ZJ'-^^ and identical topology, which 
is subject to generalized strain e and polarization stress r. Moreover, when the 
polarization stress is compatible with the difference in stiffness distribution between 
the two configurations D^'^^ and D{a) and the imposed generalized strain, the gap 
between and upper bound and the true value vanishes. 

Notice that when the optimization with respect to r is performed exactly, the 
resulting equality fl2^ is independent of the choice of the reference structure. Oth- 
erwise, e.g. when the set of admissible generalized polarizations r is constrained, 
different choices of D^^"^ generate different upper bounds on the stored energy. The 
most restrictive upper bound then corresponds to D^^^ chosen as close as possible to 
the actual generalized stiffness D{a) while maintaining the positive definiteness of 
(£>(")- £>(«))!§ To make the choice of the reference structure independent of a, it 



^Observe from Eq. (|30|) that for D'-^' — ?> -D(a), t for an arbitrary e and equality in ([29|) is 
recovered. 
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follows from the specific form of matrices -D(a) in (|22l) that the ordered eigenvalues 
of optimal Z)*-°^ are determined as the minima of the corresponding eigenvalues of 
the individual states: 



A 



(0) 



min A 



(r) 



wher e Xf' corresponds to the i-th ordered eigenvalue of matrix D^^\ cf. ( Dvorak and Srinivas , 



1999I . Section 4). 



A completely analogous procedure can be executed when selecting the reference 
structure such that (Z)*-*^-' — D{a)) becomes negative-definite, leading to a lower 
bound on the stored energy. For an indefi nite — a varia tional estimate 

of the stored energy is obtained, see e.g. iLuciano and WillisI ( 120061 ) for additional 
discussion. 



3.2. Discrete Hashin-Shtrikman variational principles 

Once inequality (12^ has been established, the Hashin-Shtrikman variational prin- 
ciples directly follow from the original energy minimization problem (124|) for a con- 
figuration a. Assume that [D^^^ — D{oi)) is positive definite and introduce a kine- 
matically admissible generalized strain e G M^'^^' obtained from a test displacement 
d G MI^^I via 

" d 



Then, we obtain the upper bound 



d 



E{d-a) 



(31) 



(32) 



+ 



- £)(a))-^' 



which yields a variational characterization of the true displacement-polarization pair 
in the form 

(33) 



''d(a), T(a)) = arg min min [/(d, r; a) 



where the Hashin-Shtrikman energy function U is defined as 

U{d, 9- a) = E^^^(d) + r^e + \t^{D^''^ - D{a)Y^' 



(34) 



with E'^^) denoting the potential energy of the reference structure (as introduced 
in Eq. (ESD). 
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Considering an arbitrary reference media and upon exchanging the order of op- 
timization, problem ( l33ll is extended into its final form: 



'd(a), rfa)) = arg stat ( min U{d,T\a 



(35) 



where the actual meaning of the "stat" operation (minimum, maximum or critical 
point) depends on the choice of the reference structure's stiffness D^^\ 

3.3. Condensed variational principle 



Similarly to IWillisI (Il977( ). we proceed with the condensation of the kinematic 
variables from the Hashin-Shtrikman function by relating the optimal displacement 
vector d to an arbitrary test polarization stress r. Due to the linearity of the problem, 
the actual value can be expressed as a superposition of two auxiliary solutions 



(36) 



where "cZ(o) denotes the polarization- independent displacement of the reference struc- 
ture subject to prescribed displacements ^d and nodal forces while is the 
displacement resulting from an internal generalized polarization stress r with = 
and = 0. The values of both components follow from the equilibrium equa- 
tions (IHl): 







"d(o) 


+ 








^d 








(0) 



(1) 



^/ - ^d. 



(37) 
(38) 



After resolving the inner optimization in Eq. 0351) . we proceed with determin- 
ing the optimal polarization. Introducing the solutions of f l37|) and f l38|) into the 
two-variable function fl35l) and exploiting the optimality conditions fl37|) and fl38|) 
yields, after some algebraic manipulations discussed in detail in Appendix \^ the 
characterization of the optimal generalized polarization stresses in the form: 



T(a) = arg stat H{t; a) 



with the condensed energy function expressed as 



(39) 



(40) 



where H^^^ = £'''°^("d(o)) corresponds to a stationary value of the potential energy of 
the reference structure, e(o) is the associated generalized strain determined from d(o) 
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and Eq. (|3T1) and r'^'^-' is the d iscrete c ounte rpart of the Green functi o n- rela ted non- 
local operators introduced by IWillid ( 119771 ) and iLuciano and WillisI (j2005[ l. linking 
the kinematic quantities to the generalized polarization stress via 

6(1) = ^"B ^ - ""B ■''B'^t = -r(°)f . (41) 

3.4- Approximation 

When considering a given probability distribution /x G A of all possible states, 
the ensemble average of the stationary energy follows from 

{Erd{a)))^ = {H{r{a)))^ £ {H{9{a)))^ , (42) 

where the shorthand notation J is used to emphasize that the actual status of the 
right hand side depends again on the choice of reference generalized material stiffness 
matrix D^^\ The equality in (H2|) remains valid since we have assumed so far that 
the optimization problem ( 139|) is resolved exactly for every configuration a. With 
the second-order description of the stochastic system at hand, a specific ansatz for 
the generalized polarization stresses given by 

2 2 



r=l r=l 



is em ployed to exploit the available second order statistics ( [T4|) optimally, cf. (IWillis . 



1973). The T^''^ and r^*^) in Eq. (gS]) are realization-independent trial and "true" 



generalized polarization stresses related to state r. 

Introducing the approximations fH5]) and the expression for the generalized ma- 
terial stiffness matrix ( 122|) into Eq. (142|) yields a variational statement in the form 



r=l 

+ ^EE-"'((x'"x'"^).-(o'°'--D")") 



r=l s=l 



Imposing the constraint /X G M(P(^^)), we can explicitly evaluate the ensemble 
averages appearing at the right hand side of Eq. (jHj). Observing that 
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as [D^'^^ —D^^^) is block-diagonal while P^^'^^ is zero-diagonal for r ^ s due Eqs. ( fTSHTBl) . 
we express the bounds/estimates in terms of the two-unit statistics 



2 



(45) 



3.5. Non-local energetic bounds and estimates 

The next step of the derivation involves the determination of the optimal state 
polarization stresses. The stationary conditions for function ( H5|) with respect to t^^^ 
variables yield a system of linear equations 

- [p(-) . - T^:^ + J2 P^"^'^ • r^^Vf^) = . e(o) (46) 

to be satisfied by the true state polarization fields r'-^'', with the subscript ~ again 
referring to the actual status of the HSW energy function. Employing the optimality 
conditions fH6|l when evaluating the terms in the square brackets in Eq. fH5l) provides 
the desired energetic bounds and estimates: 

r=l 

with the implicit dependence of r^"* on p'^^^) provided by Eq. (H6!) . The final ex- 
pression consists of the deterministic value related to the reference problem and the 
stochastic contribution, which is statistically non-local due to the incorporation of 
the discrete Green function (HT]) and P^^^^ in Eq. ( 146|) . 

3. 6. Statistics of the response 

The known values of state generalized polarization stresses allow us to obtain 
elementary statistics of additional quantities apart from the energy. For example, 
the mean value of the nodal displacements for an arbitrary /x G M(P'-^^'') follows 
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from 



(48) 



It is worth noting; t h at the prev ious relations are analogous to the results derived 
by iLuciano and WillisI (120051 12006[ ) for FE-based discretization of the HSW princi- 
ples. Therefore, adopting a similar procedure, means or conditional means of selected 
local variables can be established by post-processing the optimal state polarization 
stresses. As our focus is on the total potential ene rgy, we omit explicit expressions 
for these statistics and refer the interested reader to lLuciano and WillisI ( l2005l . 120061 ) 
for details. 



4. Examples 



Although the theory presented in the previous sections is applicable to generic 
linear discrete structures, basic features of the method are illustrated below for 
planar truss systems only. Within this context, the generalized displacement vec- 
tor and the material and kin ematic matrices introduced in Section 12.11 specialize 
to (jjirasek and Bazantl . 120011 ) : 



de 



[ 



Ul,e Vi^e M2,, 

X2,e yi,e - y2,e X2,e - 3^1,6 Z/2,e 



yi, 



] 



(49) 



where, in accordance with Figure [2], Xi^e and yi^e denote the coordinates of the e-th 
element's nodes, Uj g and Uj^e are the corresponding displacements, ig is the element 
length, Ee stands for the Young's modulus, possibly randomized in a binary sense, 
and Ae for the cross-section area. The associated generalized element strain and 
stress Se are defined as the bar's elongation and the axial force Se, respectively. 



Assuming that E^^^ > Ee^' for e 



.(2) 



1,2,..., iVg, we introduce a dimensionless 



parameter in the form 



7{±(p("))-i7(2) 



o<C<i, 



(50) 
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Figure 2: Definition of tlie geometry and kinematics of a truss element. 



to aid the visuahzation of results. The value of (± indicates the relative difference 

between the estimated mean energy of a random system described by the P^^^^ 
statistics and of deterministic structures, which solely consists of states r = 1 or 
r = 2, yielding the deterministic energies 

The regular lattice structure appearing in Figure [3] is employed as a benchmark 
problem. The prescribed boundary conditions include uniform tension and bending 
scenarios, imposed either using nodal forces (force control, FC) or by prescribed 
nodal displacements (displacement control, DC). 

The results for second-order bounds presented hereafter are accompanied by the 
statistically local counterparts of the Voigt and Reuss type, determined for a deter- 
ministic problem with spatially variable Young's moduli in the form 

= p'^^Ei'^ + il-p^^^)Ei'\ (51) 

Finally, ensemble averages obtained by direct Monte-Carlo (MC) simulations with 
= 10, 000 realizations are included reference. 

4.1. Structure with independent elements 

To demonstrate the performance of the non-local energetic bounds and estimates, 
we start with the analysis of the simplest stochastic system with a closed-form ex- 
pression for the two-unit statistics. To this end, consider a binary structure with 
the first phase assigned to each element independently with a probability 0. The 
associated second-order statistics are 

p(ii) ^ f <P if e = e', 
1 d)'^ otherwise, 
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X 




Figure 3: Sketch of the lattice structure used as a benchmark problem. 



the system is therefore statistically uniform in the sense of ( 120|) . Notice that for 
the particular network shown in Figure [3l the cardinality of the ensemble space is 
|§| = 2''^ = 4.72 ■ 10^^, whereas the computable two-unit statistical characterization 
involves ^^^^ ■ 72 = 2, 628 independent parameters only. 

The contrast of phase stiffnesses is set to E'^^^ : i?^^) = 10 : 1; the most restrictive 
upper Hashin-Shtrikman- Willis bound (HSW+) therefore corresponds to a structure 
with Ef^ = E^'^^ and the optimal lower bound HSW_ is obtained by setting Ef^ = 
E^'^\ In addition, we include variational estimates corresponding to the reference 
structure with the generalized stiffness set to the Voigt upper bound (HSW^v) or 
the the Reuss bound (HSW^r). 

Figure m presents the relative energies as a function of for all considered loading 
scenarios. The MC curves correspond to 99% confidence intervals, thereby indicating 
that a sufficient number of simulations was used to obtain reliable reference data. 
Observe that, due to the adopted dimensionless representation, the shape of the 
resulting curves is almost independent of the loading mode (tension or bending). 
For all simulations, the HSW bounds substantially narrow the domain defined by 
the first-order Voigt (V) and Reuss (R) bounds while preserving the concave/convex 
dependence of the mean potential energy on the parameter for the displacement- 
or force-driven loads, respectively. The increase in accuracy as a result of considering 
non-local spatial statistics is especially pronounced for small and large values of 0, 
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0.2 0.4 0.6 
probability ^ 

(a) tension, FC 



0.2 0.4 0.6 0.8 
probability 

(b) tension, DC 




0.2 0.4 0.6 
probability (/> 

(c) bending, FC 



0.2 0.4 0.6 0.8 
probability 

(d) bending, DC 



Figure 4: Energetics of the random truss structure with independent elements; (a)-(b) tension, 
(c)-(d) bending. The upper (HSW+) and lower (HSW_) Hashin-Shtrikman- Willis bounds are 
compared with the first-order Voigt (V) and Reuss (R) counterparts, variational estimates for the 
reference structure corresponding to the Voigt (HSW.^v) and the Reuss (HSW.^r) bound and direct 
MC simulations under displacement (DC) and force (FC) controlled loading. 
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for which even the asymptotic behavior seems to be exactly reproduced by the lower 
bound for the force-controlled conditions and the upper bound in the kinematically 
constrained case. This observation is in an agreement with analogous results for 
the Hashin-Shtrik man bounds fo r small values of the volume fractions repor ted for 
isotropic elasticity f lRoscod . Il973[ ) or scalar problems (IMilton and Nesil . Il999l ) . Even 
better agreement between the HSW predictions and MC simulations can be reached 
by employing the variational estimates generated by setting the reference stiffness 
matrix equal to the first-order bounds. Observe that for the particular system and 
loading conditions considered here, the HSW^r and HSW^v estimates effectively 
bound the reference data from above and below. As we are currently unable to 
establish optimality of the HSW+ and HSW_ bounds, further generalization of this 
fact remains an open question. 

The "energetic" displacements ( HHil as predicted by the HSW bounds appear 
in Figure [51 For both the displacement and force driven cases, the average deformed 
shapes are well reproduced by the bounds with, for example, the upper energetic 
bound corresponding to smaller displacement values. This is fully consistent with 
the fact that the upper bound represents the stiffest response for a given load and 
two-unit statistics. Moreover, the results comply with the intuitive trend suggesting 
a smaller displacement range for a larger number of kinematic constraints. 



(a) DC, maximum displacement 
scaled to 2 m. 



(b) FC, maximum displacement 
scaled to 0.25 m. 



Figure 5: Deformed structures; (a) tension, = 0.3, (b) bending, = 0.6 subject to displace- 
ment (DC) and force (FC) driven loading; displacements corresponding to lower/upper bound are 
indicated by dashed/full lines and the results of Monte-Carlo simulations are depicted in gray. 
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J^-.2. Structure with spatially variable stijfness distribution 

Having demonstrated the accuracy and hmitations of the bounds and estimates 
for a system with independent units, we focus in the second example on the effects 
of phase contrast and spatially variable distribution of material properties. For 
simplicity, the statistical model is chosen to be identical to the one introduced in 
the previous section; whereas the Young's moduli of individual elements are now 
assumed in the form 

^ l + (p-l)^, e = l,2,...,iV„ (54) 



^(2) ' 4 

where = [xi^e + X2,e)/'^ denotes the x-coordinate of the e-th element's center, 
recall Figure [3l and p > 1 is the maximum phase contrast. 

The resulting energetic bounds and estimates appear in Figure |6]for the structure 
subject to bending for low (p = 5) and high (p = 500) phase contrasts, respectively. 
We observe that in the low-contrast case, the region defined by the HSW_ and HSW+ 
bounds is substantially smaller than in the previous example, whereas the HSW^r 
and HSW^v estimates almost coincide with the MC simulation data. For increasing 
values of p, the accuracy of the bounds and estimates deteriorates; still the statis- 
tically non-local quantities provide an improvement over the first-order approaches. 
Moreover, analogously to Section 14. ![ they correctly reproduce the asymptotic be- 
havior for — and </>—)■ 1. 

^.3. Energetics of a weakest link fracture model 

In the last example, we study the ability of the bounds and estimates to cap- 
ture the energetics of a simple randomized discrete fracture model with two-unit 
probability matrices following from comp uter-generated data. 

In the deterministic setting, following iFrancfort and Marig 3 fll993[ ). a damaged 



lattice is understood as a two-state system with a "healthy" state {x^P = 1, = 0) 
and a fully-damaged state {x'e^^ = 0, x?'* = !)• As the fully damaged elements cannot 
sustain any stress, we have -E^^^ : E^^^ = 0. The transition from state 1 to state 2 
occurs when 

He = \djK,d, > Ke, (55) 

where Kf. = BjDeBe and denote the stiffness matrix and an energy threshold 
of the e-th element, respectively. 

The stochastic damage evolution is simulated using a simple combination of an 
"event-by-event" strategy and a direct Monte-Carlo algorithm as outlined in Tabled! 
For every realization of energy thresholds, a deterministic time-stepping procedure is 



20 



vj> 0.8 




0.4 0.6 
probability 

(a) FC, p = 5 



0.2 0.4 0.6 0.8 
probability 

(b) DC,p = 5 




0.4 0.6 
probability 

(c) FC,p = 500 



0.2 0.4 0.6 0.8 
probability 

(d) DC,p = 500 



Figure 6: Energetics of the random truss structure with spatially variable stiffness distribution sub- 
ject to bending; (a)-(b) phase contrast p = 5, (c)-(d) phase contrast p = 500. The upper (HSW+) 
and lower (HSW_) Hashin-Shtrikman- Willis bounds are compared with the first-order Voigt (V) 
and Reuss (R) counterparts, variational estimates for the reference structure corresponding to 
the Voigt (HSW^v) and the Reuss (HSW^r) bound and direct MC simulations under displace- 
ment (DC) and force (FC) controlled loading. 
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1. Input: probability distribution of thresholds Ke for e = 1, 2, . . . , A'^e, loading 
scenario f{t)d with /(O) = 0, /(I) = 1 and / being monotonic, number of 
realizations and sampling times = Tq < ri < . . . < rj\/ = 1 

2. For z = 1, 2, . . . , iV 

(a) initialize = 0, r = and x^^in) = 1 G for A; = 1, 2, . . . M 

(b) set Ee = E and generate thresholds Kg for e = 1, 2, . . . , A^'e 

(c) Repeat 

i. set '^'X^^'^^Tk) = and = for e E (f) and j E t 

ii. for current values of E^, compute element energies (15 5 p using dis- 
placements dg determined from (|H]) with "^d = d; e = 1,2, . . . , 

iii. set £ = {e = 1, 2, . . . , : > 0} 

iv. compute = /H^ for e G £ 

V. set t = minege{te} 

vi. set (p = {e E s : te = t} 

vii. set T = {k = 1,2, . . . , M : Tk > t} 
until t > 1 or £ = 

3. using samples {x^'K'^k)x^''^'^ iTk)}f=i, compute P^lcVfc) ior k = 1,2, . . . , M 

4. using samples {ifW(rfc)}^^, compute 'H^*"(rfc) for k = 1,2, . . . , M, where 
H^'Hn) = Ef=i^e, with He provided by using displacements d^ deter- 
mined from ([8]) with ^d = f{Tk)d and E^ = Xe\Tk)E 

Table 1: Conceptual implementation of the weakest link algorithm. 
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executed. First, the structure is subjected to the full displacements and the energies 
corresponding to individual elements are extracted (step ii.). Then, in steps iii.-vii., 
the element (s) with the highest damage affinity is (are) determined and removed 
from the structure (step i.). This is accompanied by the bookkeeping of element 
deactivation times in terms of the auxiliary matrices x'*^(Tfc), k = 1, 2, ... , M. The 
procedure is repeated until the imposed loading sequence is completed or until the 
structure loses its integrity; cf. the termination conditions of step (c). After the 
sampling phase is complete, quantities of interest related to the damage statistics 
follo w from post-processing of the simulation results (s t eps 3. and 4 . ); see also 



e.g. ( ISharif-Khodaei and Zemanl . l2008l : lAlava et al.l . l2006l : IVellinga et al.l . 120081 ) and 



references therein for more details. 

In particular, we consider the lattice shown in Figure |3] subject to an imposed 
displacement, which is parametrized as 

u{t) = u^^^Vi, < t < 1, (56) 



leading to a linear scaling of the element energies fl55|) with respect to the pseudo- 
time t. 

For simplicity, identical values of the Young's moduli, E^^^ = E, are assumed 
in the healthy state for every element and the thresholds are taken as uniformly 
distributed independent random variables ranging from k~ = ^EeAe{A£'^J)'^ /if. to 
with the critical elongation A£™ = 10 ^le- The target displacement is set 
to Mmax = 10~^ for uuiform tension, whereas for the bending case Wmax = 1-5 • 10~^. 
The time evolution of the system was sampled at M = 100 uniformly distributed 
time instants using = 10, 000 realizations. 

Snapshots of the damage progress quantified in terms of one- and two-unit proba- 
bilities appear in Figure [7] for both the uniform tension and bending scenarios. Note 
that the diagonal entries of the P^^^\t) matrix provide, in the current setting, the 
survival probability for each element at time t. The off-diagonal values correspond 
to the simultaneous survival probabilities for two elements. 

In the case of uniform tension, the damage distribution appears to be rather 
diffuse during the whole loading procedure. Low survival probabilities are localized 
in the horizontal links in the corners of the network and in the structure's interior, as 
these exhibit the highest failure affinity for the deterministic solution with 2 

(pseudo-time t = 0.1) |f| The slight non-symmetry of the one-unit probabilities 
visible in Figure [7] results from the finite number of realizations used to sample 



^Therefore, in the terminology introduced by Eq. ()20[) . the initially statistically uniform structure 
evolves to a statistically non-uniform system. 
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tension bending 

f =0.1 





/ =0.6 




20 30 40 50 

element number 



20 30 40 50 
element number 



Figure 7: Evolution of the one- and two-unit probability matrix entries as determined by 

Monte-Carlo simulations. The snapshots for t = 0.1 and t = 0.6 show one-unit probabilities, 
whereas complete two-unit probability matrix i^^lotted for t = 0.9. 



the statistics. With increasing load, the remaining horizontal elements start to fail 
due to the energy redistribution after failure of the first group of weakest links, 
leading to a decrease of the survival probabilities for the interior elements {t = 0.6). 
In the case of bending, on the other hand, the damage is triggered mainly in the 
upper and lower rows of elements and subsequently propagates to internal links. 
Additional details on the simultaneous survival probabilities are given by the full 
two-unit probability matrices, shown in Figure [7] for t = 0.9. Recalling element 
numbering introduced in Figure [3l the plots reveal, among other things, simultaneous 
survival probabilities close to one for vertical element groups 21-24 and 37-40, which 
connects the nodes where the load is imposed, and horizontal elements 9-10, at the 
neutral axis, in the case of bending. It follows from these results that the second-order 
statistics correctly capture the dominant fracture mechanisms, thereby providing a 
well-founded statistically non-local damage parameter for stochastic damage theories. 

The evolution of the normalized stored energy ( is depicted in Figure El The 
choice of reference structures to generate the bounds or estimates is identical to that 
in the previous example. Note that the lower bounds and the HSW^r estimate are 
identically equal to zero due to the zero stiffness assigned to the weaker state. The 
overall character of the energy evolution is consistent with the previous discussion. In 
particular, two failure modes are clearly visible for the bending problem, whereas the 
gradual decrease in energy for uniform tension corresponds to more a diffuse damage 
character. Similarly to Section 14. H the first-order Voigt bounds correctly predict the 
trend of the energy reduction, but under predicts its magnitude. The HSW bound 
remains highly accurate until t ^ 0.2. A certain discrepancy, however, is observed for 
the limit value of the dimensionless energy, e.g. = 0.24 instead of the reference 

^Mc ^ g uniform tension, cf. Figure [8]|^a). This difference again arises due to the 
fact that the upper HSW bound corresponds, in the ensemble average sense, to 
the response of the stiffest structure compatible with the two-unit statistics shown 
in Figure [71 for which the minimum probability reaches ~ 30% for both loading 
sequences. The results of the MC simulations, on the other hand, represent one 
particular stochastic system implicitly defined by the randomization procedure and 
the event-by-event solution algorithm. Nevertheless, the added value of the second- 
order relations becomes apparent when considering the Voigt prediction = 0.74. 
Slightly more accurate values of the residual relative energy can be generated using 
the HSW^v estimate, for which we obtain (^^sw^v = 0.21. The general character of 
these conclusions is further supported by the analogous results of the bending mode 
shown in Figure ElJ^b). 
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Figure 8: Energetics of the weakest link model; (a) tension, (b) bending. 



5. Conclusions 

In this work, variational bounds and estimates of the Hashin-Shtrikman- Willis (HSW) 
type for two-phase random structures have been derived and verified against the re- 
sults of direct Monte-Carlo (MC) simulations. The most important findings can be 
summarized as follows: 

• When applied to discrete structures instead of continua, the derivation of the 
HSW principles becomes reasonably straightforward and requires only elements 
of matrix structural analysis and linear algebra. 

• The variational framework naturally incorporates general statistically non- 
uniform systems, including the specification of the set of all admissible two-unit 
probability matrices. Thus, the developed statistically non-local theory is able 
to treat the statistics of systems with independent and highly correlated units 
in a unified manner. 

• The second-order bounds and estimates provide a computationally feasible al- 
ternative to direct MC simulations. 

• When applied to the displacement-driven damage problems, the upper bound 
on the stored energy delivers a reasonably accurate approximation to the total 
potential energy at the damage initiation and propagation without any ad- 
justable parameters. When complemented with an accurate representation of 
the irreversible energy dissipation during the damage process, the scheme has 
the potential to provide a consistent variational model for damage evolution in 
discrete media. 
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The developed statistically non-local discrete theory can also provide a convenient 
starting point for the coarsening towards an equivalent continuum representation. 
The already announced applications, i.e. a variational approach to deterministic and 
stochastic damage mechanics of finite-sized lattices and incorporation of higher-order 
statistics, will explored and reported separately in future publications. 



A. Derivation of condensed energy 



Recall that the expression for the condensed energy follows from the value of the 
Hashin-Shtrikman function fl34|) with the displacement vector optimally adjusted to 
a trial generalized polarization stress r. In particular, evaluating the energy for the 
displacement determined from ( l36i) yields, in a partitioned format: 



(0) 



d^^^y^f + r^ -^B -B ] 



'd 



(0) 



"d 



A suitable re-arrangement of terms appearing in the previous relation together 
with the optimality conditions (1371) and (l38l) leads to 



by definition 



f/("d,T;a) 



''d(o) 


T r 







"(0) 
^d 



d(0) J 



=6(0) by definition 



+ f ^ [ ""B ''B ] 

=-r(0)r by Eq. glj 



d(0) 
^d 



=0 by Eq. (|38} 



^ 2 ' 



"d(i) +|"d(i)T (K(°)"d(i) + ^"B'^r) 



=0 by Eq. l(37)l 



^d 
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which coincides with the expression for the condensed energy as appearing in E g. (I40l). 
Note that similar development for the continuous case is available in, e.g., IWillisI 
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( 1977 ) and in Iwillisi ( 1981 ). with the FE-discretization treated in 

toom . 
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